Ferroelectricity and isotope effects in hydrogen-bonded KDP crystals 
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Based on an accurate first principles description of the energetics in H-bonded KDP, we conduct 
a first study of nuclear quantum effects and of the changes brought about by deuteration. Cluster 
tunneling involving also heavy ions is allowed, the main effect of deuteration being a depletion 
of the proton probability density at the O-H-O bridge center, which in turn weakens its proton- 
mediated covalent bonding. The ensuing lattice expansion couples selfconsistently with the proton 
off-centering, thus explaining both the giant isotope effect, and its close connection with geometrical 
effects. 
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Potassium dihydrogen phosphate (KH2PO4, or KDP) 
belongs to a family of ferroelectric (FE) crystals where 
molecular units are linked by hydrogen bonds, the ferro- 
electricity being connected to proton off-center ordering 
in the bonds. A characteristic feature of this family is the 
large increase in the Curie temperature T c upon deuter- 
ation. In this particular case, it goes from ~ 122 K in 
KDP to ~ 229 K in the deuterated compound (DKDP). 
The origin of this huge isotope effect is still controversial, 
and has been mostly understood in terms of the quan- 
tum tunneling model proposed by Blinc, jlj later modi- 
fied by inclusion of the coupling between proton motion 
and the K-PO4 dynamics. [2| While direct experimental 
indications of tunneling have recently emerged S, the 
connection between proton tunneling and isotope effect 
remains unclear. There is in fact strong evidence, 
that isotope substitution acts rather through a geomet- 
rical modification of the hydrogen bonds, Q with an ex- 
pansion of the O-H-0 distance. The proton off-centering, 
and thus the corresponding increase of lattice parameter 
upon deuteration, appear to be remarkably correlated to 
the increase of order parameter and of T c . These find- 
ings stimulated new theoretical work where some of these 
facts could be addressed without invoking tunneling, fij] 
however so far only at a rather phenomenological level. 

In the first part of this letter we investigate, using elec- 
tronic structure calculations within Density Functional 
Theory (DFT), the relationship between proton order- 
ing, polarization, and geometry in KDP. In the second 
part, we introduce a study of energy and subsequently the 
quantization of the collective ion displacements in small 
KDP clusters, embedded in a host paraelectric lattice. 
These calculations demonstrate the difference between 
deuterated and protonated KDP, the more delocalized 
proton bridging the oxygens and pulling them together 
more effectively than the deuteron. This phenomenon, 



which is at the root of the geometric effect, is further il- 
lustrated by solving in the last part a selfconsistent non- 
linear model. 

For the DFT calculations we use two different ap- 
proaches: one employing a basis set of confined pseu- 
doatomic orbitals (SIESTA), || another a plane wave 
(PPW) representation. || For the first we choose a 
double-zeta basis set with polarization functions, and 
an orbital confinement energy E c — 50 meV. In the 
second, we set the energy cutoff to 150 Ry. In both 
cases exchange-correlation terms are computed using 
a gradient-corrected (GGA) functional, |10| and norm- 
conserving pseudopotentials |pd| are employed to repre- 
sent the interaction between ionic cores and valence elec- 
trons. We also include the nonlinear core corrections for 
a proper description of the K ion. The Brillouin zone 
T-point alone provides a sufficient sampling in the large 
supercells used. 

The paraelectric (PE) structure of KDP is body- 
centered tetragonal with two KH2PO4 units per lattice 
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FIG. 1. Schematic view of the internal structure of KDP 
along the tetragonal axis. P and K coordinates relative to 
c along this axis are indicated. Covalent and H-bonded hy- 
drogens in the FE phase are attached to the corresponding 
oxygen through full and broken lines, respectively. 
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site. We use the conventional bet cell doubled along the 
tetragonal c axis (64 atoms) to describe homogeneous 
distortions, and a conventional fct cell doubled along 
the c-axis (128 atoms) for local distortions. The internal 
structure is depicted in Fig. 1. Above T c the protons 
occupy with equal probability two equivalent off-center 
positions in the H-bond separated by a distance S, E^] 
while below T c they order in such a way that each PO4 
group has two covalently bonded and two H-bonded hy- 
drogen atoms. 

We first analyse the relationship between proton order- 
ing and polarization. We start from the average experi- 
mental structure of the PE phase of KDP at Tf DP +5 K, 
fl2"f with the hydrogens centered in the O-H-0 bonds. By 
fully relaxing the atomic positions in the bet cell, the H 
atoms move collectively off-center towards the 02 oxy- 
gens, and away from 01, as indicated in Fig. 1. Analysis 
of Mulliken populations and charge density differential 
maps |^3) indicates that the off-centering of the H atoms 
induces a covalent charge displacement from 02 towards 
the 02-H bonds, whereas the Ol-H bonds weakens into a 
hydrogen bond. In addition, there is a charge flow from 
the 02-P to the Ol-P bonds, accompanied by an increase 
of the P-02 and a decrease of the P-01 distances. P 
atoms are thus driven off-center in the PO4 tetrahedra, 
which polarize further. Unbalanced electrostatic forces 
induce a displacement of the K + ions along the c-axis, 
towards the charge-excess (01) side of PO4 units. A de- 
tailed description of the classical FE distortion in KDP 
can be found in Ref. jl3| . These results are also in agree- 
ment with Ref. p"4| ] . In order to identify the driving mech- 
anism of the FE instability, we investigate the ab initio 
potential energy surface (PES) as a function of the pro- 
ton off-centering parameter S = doo ~ 2doH , and of the 
K-P relative displacement along c, 7 = c — 2dnp, which 
provides a measure of the polarization. |]l5|| By fully re- 
laxing the oxygen positions for each chosen (5, 7) pair, we 
obtain a two-dimensional double-well PES with a saddle 
point at S = 7 = 0, whose contours are reported in the 
inset to Fig. 2. According to this PES, the crystal is sta- 
ble against polarization (7 ^ 0) unless the protons are 
ordered off-center (<5 / 0). In Fig. 2 we show cuts of 
the PES at different values of 7, indicating that, even for 
vanishing 7, the energy minimum corresponds to a finite 
5, i.e. protons are always collectively unstable at the H- 
bond centers. Therefore, we confirm that the source of 
the FE instability is indeed the H off-centering. 

In order to address next the quantum effects, in partic- 
ular the isotope substitution, we need to introduce quan- 
tum fluctuations in the above classical picture. Barring 
for the time being a full brute force quantum mechani- 
cal calculation for all the ionic degrees of freedom of the 
infinite system, we take a different approach to the prob- 
lem. Important quantum effects are identified as those 
involving correlated H motions (as shown in Fig. 1) with 
relaxation of K and P ions, most favorable for exhibiting 
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FIG. 2. Energy profiles as a function of S, for values of 
7 = (solid line), 0.02 (short-dashed), 0.05 (dot-dashed) and 
0.1 (long-dashed) A. The minima are always at S 0, and 
for 7 > 0.02 A the curves exhibit a single minimum. The 
inset shows equispaced energy contours (step = 13.6 meV/ 
KH2PO4 unit). The minima at (5, 7) ~ ±(0.3, 0.15) °A lie ~ 50 
meV below the saddle point at (0,0). 

FE instabilities. We shall consider this correlated pattern 
as a single collective coordinate <5 C . Quantization of this 
coordinate will moreover be carried out within a finite- 
size cluster. Even though quantum fluctuations span in 
principle all length scales, those that occur at short range 
should be sufficiently revealing at least away from critical 
points. 

In this spirit, we consider a series of small KDP clusters 
of increasing size, embedded inside an undistorted host 
PE matrix. For these clusters, we first determine, clas- 
sically, the total energy variation as a function of S c /2. 
The clusters comprise N hydrogens, in the following or- 
der : (a) N=l; (b) N=4, fully connecting a PO4 group 
to the host; (c) N=7, connecting two PO4 groups; (d) 
N=10, connecting three PO4 groups. In order to ascer- 
tain the effect of the volume increase seen upon deuter- 
ation, we repeat the calculations by expanding the host 
structural parameters to the corresponding experimental 
values of DKDP at Tf KDP +b K. |[2| The results show 
(Fig. 3, solid curves) that there is an instability only 
when the cluster size is sufficiently large, thus providing 
a measure of the FE correlation length. The instabil- 
ity is much stronger (note the larger energy scale), and 
the correlation length accordingly shorter, with the ex- 
panded structural parameters of DKDP. This is in close 
agreement with the experimental trend, showing that the 
FE order grows with volume. Finally, the involvement of 
P and K motions is important. The dashed curves in 
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Fig. 3 represent the energy of a distortion where only H 
atoms are allowed to move, all other atoms being kept 
fixed. The much higher energies indicate that, not sur- 
prisingly this kind of oversimplified distortion is far from 
real, leaving, e.g., the KDP lattice fully stable. 

The next step in order to study the quantum effects 
due to proton or deuteron zero point motion, is to quan- 
tize the clusters with respect to the local collective coor- 
dinate 6 c /2. The corresponding kinetic energy involves 
the mass of each ion proportionally to the square of its 
displacement, and will change upon deuteration. We find 
that the deuterium (hydrogen) effective mass for this cor- 
related motion is about /i£> = 3.0 (/ijj = 2.3) proton 
masses in DKDP (KDP) clusters. Solving Schrodinger's 
equation for the single variable S c /2 with mass hh,d and 
effective potentials of Fig. 3 is trivial. The ground state 
(GS) energy levels which lie below the energy barrier at 
<5 c /2 = are indicated by dotted lines. A negative en- 
ergy signals the occurrence of tunneling between + and 
- 8 c /2, and its onset provides a rough indication of the 
correlation volume. The results show that it comprises 
more than N=10 hydrogens in KDP, but no more than 
N=4 deuteriums in DKDP. In both systems, we eventu- 
ally expect the collective tunnel splitting to converge to 
zero as N — * oo, signaling long-range ferroelectricity; the 
tunnel splitting should conversely remain nonzero in the 
quantum paraelectric state, attainable at high pressure. 
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Remarkably, deuteration at fixed geometry does not 
affect the fact that the tunnel splittings in large clus- 
ters [ [l6| arc exceedingly smaller than ksT c , which ren- 
ders T c practically insensitive. In addition, the effective 
mass change is reduced by the involvement of the heavy 
nuclei. This agrees with high-pressure experiments, PJl7| 
where under fixed structural conditions the isotope effect 
appears to be only very slight. 

The fact that energy barriers in DKDP are much larger 
than those in KDP implies that quantum effects are sig- 
nificantly reduced in the expanded DKDP lattice. Af- 
ter this observation, the reason for the geometric effect 
becomes clear: in all KDP clusters considered here the 
wavefunction (WF) for the collective coordinate 8 C is 
much more delocalized, with a larger amplitude near the 
H-centred position S c = between 01 and 02, than in 
DKDP. In KDP in other words, zero-point motion pushes 
protons towards the center, much more effectively so than 
deuterons in DKDP. An increased probability for the pro- 
ton to bridge midway between the two oxygens is not ir- 
relevant to the crystal cohesion, and we identify precisely 
that as the element which compresses the cell from a 
larger classical value to the smaller value found for KDP. 
To estimate an upper limit to that effect, we compare 
the lattice parameter and the O-H-0 bridge length doo 
of two classical electronic calculations: one with the hy- 
drogens forced to sit in the bridge center, the other with 
H fully off-center in the classical FE state. The result is 
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FIG. 3. Ab initio energy profiles for classical distortions of 
clusters embedded in undistorted PE structures of (a) KDP 
and (b) DKDP. Reported are clusters of: N=4 (diamonds), 
N=7 (squares), and N=10 (circles). Full symbols and solid 
lines: P and K positions inside the cluster are allowed to 
relax as H are off-centered; empty symbols and dashed lines: 
only H atoms are moved. Dotted lines: GS energies (only 
negative ones, signaling tunneling). 

4.42 A when H is centered. At the equilibrium volume, 
we estimate that proton centering creates an equivalent 
pressure of approximately 20 kbar. Thus, a centered hy- 
drogen acts as a very strong attraction center that pulls 
the two bridge oxygens together, the whole lattice vol- 
ume shrinking by 2.3 % as a result. In the true high- 
temperature PE phase the hydrogens are of course not 
exactly centered, but appear on both sides of the H-bond, 
thus reducing the magnitude of the effect. 

Deuteration of KDP to DKDP at fixed geometry re- 
sults in a depleted probability for the deuterium to reach 
the bridge center, as shown in Fig 4(a) where we plot 
the proton and deuteron WF in the DKDP potential of 
the N=7 cluster of Fig. 3(b). This depletion releases a 
fraction of the O-H-0 bond grip, causing a small lattice 
expansion which has the effect of increasing the poten- 
tial wells depth, making the deuteron even more local- 
ized, and so on in a self-consistent manner. The overall 
self-consistent effect is eventually much larger than the 
mere fixed-geometry replacement of the proton with the 
deuteron mass. In fact, for the same cluster size in KDP, 
instead of a double-peak, the WF exhibits a broad single 
peak centered in 8 c /2 = (see Fig. 4(a)). This mecha- 
nism is now capable of explaining, at least qualitatively, 
the increase in the order parameter and T c , which are re- 
lated to the mass only indirectly, through the geometric 
modification of the length and energy scales. 

We constructed a simple model that demonstrates the 
nonlinear behaviour arising from isotope substitution in 
KDP by adding to the effective hydrogen potential in the 
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FIG. 4. WF in the N=7 cluster PES for (a) ab initio and 

(b) self-consistent model calculations. Solid (dashed) lines are 
for DKDP (KDP). Dotted line is for H in the DKDP PES. 

(c) WF peak separation S p as a function of the effective mass 
fi (given in units of the proton mass) for the self-consistent 
model (circles) and for fixed DKDP potential (squares). Lines 
are guides to the eye. 

cluster Schrodinger equation a quadratic WF-dependent 
term of the form V e s(x) = Vo(a;) — fc| 1 3/(a;) | 2 , where x = 
S c /2 and Vo(x) is a quartic double- well similar to those of 
Fig. 3. The nonlinear ^(a;)) 2 term signifies the positive 
feedback discussed above, since a decreasing mass will 
increase ^(x) | 2 at the center, which will in turn lower the 
barrier at that point, further increasing ^(x)) 2 , and so 
on. By choosing appropriately k and the "bare" potential 
Vo (x) , we find that the corresponding ab initio profiles for 
KDP and DKDP can be qualitatively reproduced (Fig. 4 
(b)). The self-consistent solution evolves from a double- 
peak to a single-peak situation when changing the mass 
from "pure DKDP" to "pure KDP" , as shown by the 
circles in Fig. 4 (c). Such a large mass dependence, in 
striking contrast with the very weak dependence obtained 
at fixed DKDP potential and geometry (squares), can 
now explain the large and controversial isotope effect in 
the ferroelectricity of KDP. 

Summarizing, we confirm that the H off-center order- 
ing is the source of FE instability. First principles calcu- 
lations in a model PE phase show that distortions involv- 
ing only H atoms do not display significant instabilities. 
Instead, heavy cluster tunneling involving also K and P 
ions, explains the "double occupancy" in the PE phase, in 
agreement with the observation of P-ion double-site dis- 
tribution. |D| Although tunneling is of course the crucial 
ingredient in the isotope effect, its mere change at fixed 
lattice geometry accounts only for a very small fraction 
of the full isotope effect. The main effect of replacing 



deuterons with protons appears to be an enhancement of 
the quantum delocalization of the proton in the bond cen- 
ter region. That gives rise to a strong proton-mediated 
covalency in the O-H-0 bridge, which shrinks the lat- 
tice, which further delocalizes the proton, and so on in 
a nonlinear loop. In the end, the huge isotope effect on 
T c is dominated by the geometrical effect. This effect, 
however, is triggered by tunneling, thus reconciling these 
two aspects which were largely debated in the past. 
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